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Stability of toroidal magnetic field in a stellar radiation zone is considered for the cases of uniform and differential rotation. 
In the rigidly rotating radiative core shortly below the tachocline, the critical magnetic field for instability is about 600 G. 
The unstable disturbances for slightly supercritical fields have short radial scales ~1 Mm. Radial mixing produced by the 
instability is estimated to conclude that the internal field of the sun can exceed the critical value of 600 G only marginally. 
Otherwise, the mixing is too strong and not compatible with the observed lithium abundance. Analysis of joint instability 
of differential rotation and toroidal field leads to the conclusion that axisymmetric models of the laminar solar tachocline 
are stable to nonaxisymmetric disturbances. The question of whether sun-like stars can posses tachoclines is addressed 
with positive answer for stars with rotation periods shorter than about two months. 
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1 Introduction 

The problem of stability of toroidal magnetic fields in stel- 
lar interiors first addressed by Tayler (1973) is highly rele- 
vant for the dynamics of stellar radiation zones. The stabil- 
ity controls which primordial fields can survive inside the 
stars on evolutionary time scales. Magnetic instabilities are 
important for transport of angular momentum and chemi- 
cal species in the radiation zones (Barnes, Charbonneau & 
MacGregor 1999). The stability problem is also relevant to 
the origin of the solar tachocline (Gilman 2005). 

This paper concerns the stability of toroidal magnetic 
fields in rotating radiation zones. The analysis is global in 
horizontal dimensions but the radial scales of unstable dis- 
turbances are assumed as short. The smallness of the ra- 
dial scales is a consequence of the stable stratification. Our 
mathematical formulation is close to the thin layer approx- 
imation of Cally (2003) and Miesch & Gilman (2004). Our 
aims, however, are quite different. We do not address the 
growth rates of highly supercritical disturbances but focus 
on the marginal stability. We also include finite diffusion 
that is necessary when stability of rigidly rotating region 
beneath the tachocline is concerned. When diffusion is ne- 
glected, the most rapidly growing disturbances have indefi- 
nitely short radial scales and the instability can be excited by 
indefinitely small magnetic field. The unphysical zero wave 
lengths and field amplitudes are prevented by finite diffu- 
sion. The marginal field strength for the rigidly rotating re- 
gion just beneath the tachocline is about 600 G. Estimation 
of the radial mixing produced by the instability suggests that 
the marginal value is close to the upper bound for the inter- 
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nal field. The mixing is too strong and is not compatible 
with the observed lithium abundance if the marginal value 
is considerably exceeded. 

Differential rotation can drive its own instability even 
without magnetic fields (Watson 198 1 ; Dziembowski & Ko- 
sovichev 1987). The instability has been extensively studied 
in 2D approximation of strictly horizontal disturbances. An 
account for radial displacements in our 3D formulation re- 
duces the magnitude of differential rotation required for the 
hydrodynamic instability. 

We find that axisymmetric models of a laminar solar 
tachocline do not suffer from nonaxisymmetric magnetic in- 
stabilities. The question arises whether stars other than the 
Sun can possess tachoclines. Evidences on stellar differen- 
tial rotation used in combination with laminar tachocline 
theory provide a positive answer for all solar-like stars 
whose rotation is not too slow. 

2 The model 

The background state of our model includes differential ro- 
tation and toroidal magnetic field. The angular velocity dis- 
tribution is parameterized as 

O = n Q (1 -a cos 2 9) , (1) 
where flo is the equatorial value of ft and 9 is colatitude. 
The toroidal field is expressed in terms of the Alfven angu- 
lar frequency fi A , i.e. 

B = r sin 9^/JIop Oa (r, 9) e$, (2) 
where is the azimuthal unit vector. The field is assumed 
to consist of two broad belts antisymmetric about the solar 
equator: 

SI A = flob cos 9. (3) 
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Fig. 1 (online colour at: www.an-journal.org) The buoyancy fre- 
quency (Eq. 4) in the upper part of the radiative core of the Sun 
after the model by Stix & Skaley (1990). The convection zone of 
the model includes the overshoot layer so that jV/fi© is large im- 
mediately beneath the convection zone. 

b is the normalized field strength, 6=1 means approximate 
equipartition of rotational and magnetic energies. 

The equations for linear perturbations written in the Ap- 
pendix were derived in Kitchatinov & Riidiger (2007). The 
equations apply to the disturbances that are global in hor- 
izontal dimensions but short-scaled in radius. The reason 
why the radial scales are short is the stable (subadiabatic) 
stratification. The buoyancy frequency, N, 
g ds 



O.OO 0.02 0.04- 0.06 



N 



C p 5r' 



enters the normalized equations via the parameter 



A = 



N 
Qokr ' 



(4) 



(5) 



where g is the gravity, C p is the heat capacity at constant 
pressure, s = C v In (P/p 7 ) is the specific entropy, r is the 
(heliocentric) radius and k is the radial wave number of the 
disturbances. The ratio N/CIq in stellar radiation zones is 
a large quantity (Fig. 1). After Eq. (5), disturbances can 
avoid the stabilizing effect of stratification by decreasing 
their radial scales. Our results shall confirm this short-wave 
approximation with kr <C 1. 

The equations for small disturbances listed in the Ap- 
pendix allow two types of equatorial symmetry. We use the 
standard notations Sm and Am for symmetric and antisym- 
metric modes, to is the azimuthal wave number. Sm modes 
have mirror-symmetry (about the equatorial plane) distur- 
bances of the magnetic field B' (symmetric B'^, B' r and 
antisymmetric B' g ), mirror-antisymmetric velocity u' (sym- 
metric u' g and antisymmetric and u' r ) and antisymmet- 
ric perturbations of entropy s' . The symmetries reverse for 
Am modes. These symmetry conventions are the same as in 
Gilman, Dikpati & Miesch (2007). 

3 Results and discussion 

3.1 Beneath the tachocline (a = 0) 

Helioseismology shows that the rotation below the 
tachocline is almost rigid (Schou et al. 1998). The uniform 



Fig. 2 (online colour at: www.an-journal.org) Stability map for 
ideal fluid of zero diffusion and rigid rotation. Indefinitely small 
radial scales are preferred. 



rotation has minimum kinetic energy for given angular mo- 
mentum. Rotational energy cannot, therefore, feed any in- 
stability so that magnetic instabilities can only be consid- 
ered. 

3.1.1 Ideal fluids 

Microscopic diffusion in radiative cores is small and often 
neglected. We start with the case of ideal fluids with diffu- 
sion parameters equal zero. The stability map for this case 
is shown in Fig. 2. 

Only the kink-type (to = 1) disturbances are found un- 
stable. The instability remains active for indefinitely small 
field amplitude b though the range of radial scales of un- 
stable excitations reduces with decreasing b. Therefore, ro- 
tation does here not suppress the Tayler instability. The 
character of rotational influence crucially depends on the 
toroidal field profile. Pitts & Tayler (1985) found strong ro- 
tational suppression but suggested that it might be restricted 
to the case in which the ratio of rotational to Alfven veloc- 
ities is uniform. The results of Fig. 2 confirm the expecta- 
tion. The results were found with same numerical code as 
formerly reported strong rotational quenching for uniform 
fl A (Kitchatinov & Riidiger 2007). 

The ideal instability for very weak fields is, however, of 
mainly academic interest because the unstable disturbances 
have indefinitely small radial scales. Finite diffusion must, 
therefore, be included. 

3.1.2 Finite diffusion 



From now on the finite values 



2 • 10 



-10 



4-10" 



10' 



(6) 



that are characteristic of the upper part of the solar radiative 
core are used for the diffusion parameters (Eq. A6). The 
relation x 3" V ^ v of Eq. (6) is typical of stellar radiation 
zones. 

The stability map computed with finite diffusion is 
shown in Fig. 3. Small radial scales are now stable due to 
finite magnetic diffusivity. 
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Fig. 3 (online colour at: www.an-journal.org) Stability map for 
finite diffusion and rigid rotation. 
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Fig. 4 (online colour at: www.an-journal.org) Growth rates in 
units of Qq for b = 0.01 (left) and 6 = 0.1 (right). The largest rates 
belong to A 5; 0.1. All curves for most rapidly growing modes. 

For the upper radiation zone of the Sun with p ~ 
0.2 g/cm 3 it follows 

~ 10 5 6 G, A ~ 10A Mm. (7) 
After Fig. 3, the critical magnetic field for the onset of the 
instability is thus about 600 G. The modes that first become 
unstable when this field strength is exceeded have vertical 
wavelengths between 1 and 2 Mm. For larger field strengths, 
there is a range of unstable wavelengths. The maximum 
growth rates remain, however, at the wavelengths ~ 1 Mm 
(Fig. 4). The growth rates are rather small. 

The instability radially mixes chemicals. It can thus be 
relevant for the radial transport of light elements (Barnes 
et al. 1999). The effective diffusivity, Dt ~ A 2 a can be 
estimated from our linear computations. With Eq. (7) this 
yields 

D T ~ 10 12 A 2 a cm 2 /s, (8) 
where a is the normalized growth rate given in Fig. 5. Diffu- 
sivities in excess of 10 4 cm 2 s _1 in the upper radiative core 
are not compatible with the observed solar lithium abun- 
dance. Hence, the toroidal field amplitude can only slightly 
exceed the marginal value of about 600 G. 

The estimation of light elements mixing seem to exclude 
the possibility of hydromagnetic dynamo driven by mag- 
netic instability in the solar radiation zone (Spruit 2002). 
This possibility is also not supported by direct numerical 
simulations of Zahn, Brun & Mathis (2007). 
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Fig. 5 (online colour at: www.an-journal.org) The growth rate 
a in units of fio as function of the toroidal field amplitude. The 
dotted line shows the parabolic approximation a = 0.1& 2 . The 
scale on the right gives the radial diffusivity of chemical species 
estimated after Eq. (8). 

3.2 Inside the tachocline (a > 0) 

The helioseismology detections of the tachocline central ra- 
dius r c = (0.693 ± 0.002) R and (equatorial) thickness 
Ar = (0.039±0.013)R S (Charbonneau et al. 1999b) place 
the tachocline mainly if not entirely beneath the base of 
convection zone at r = 0.713 Rq (Basu & Antia 1997; 
Christensen-Dalsgaard, Gough & Thompson 1991). The 
stability analysis of toroidal fields of radiative core can thus 
be applied to the tachocline. The differential rotation with 
a > in Eq. (1) must also be included. 

3.2.1 Hydrodynamic stability 

Differential rotation can be unstable even without magnetic 
fields. The hydrodynamic stability of the tachocline has 
been studied in 2D approximation of strictly horizontal dis- 
placements (Charbonneau, Dikpati & Gilman 1999a; Ga- 
raud 2001). Here we extend the treatment to 3D. 
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Fig. 6 (online colour at: www.an-journal.org) The map of hy- 
drodynamic instability of differential rotation (1). Most unstable 
are the perturbations of Al symmetry type with the vertical scale 
A ~ 0.6. The critical magnitude of latitudinal shear is reduced to 
0.21 compared to the 0.28 of 2D theory. 
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Fig. 7 Isolines of small growth rates on the plane of differential 
rotation and toroidal field parameters a and b of Eqs. (1) and (3). 
Laminar tachocline models are stable. 

The resulting stability map is shown in Fig. 6. Al and S2 
modes can be unstable for positive a. The marginal values 
of a approach constants for large A. The large- A limit repro- 
duces the results of the 2D theory as it should (Kitchatinov 
& Riidiger 2007). The minimum a- value for the instabil- 
ity corresponds, however, to finite vertical scales of about 6 
Mm. 

3.2.2 Joint instabilities 

Growth rates of joint instabilities for b ~ 1 have been dis- 
cussed by Cally (2003) and Gilman et al. (2007). Now we 
concern the marginal stability that is important in relation to 
the models of laminar tachocline formed by a weak internal 
magnetic field (Riidiger & Kitchatinov 2007). 

Isolines of the very small growth rates for Al and SI 
modes for the a and b parameters of differential rotation and 
magnetic field are shown in Figs. 7 and 8. The growth rates 
are optimized in A, i.e. A has been varied to find the max- 
imum growth rates for given a and b. From Sect. 3.2.1 we 
know that differential rotation is stable to SI disturbances. 
For sufficiently large a, the S 1 modes become unstable for 
very small magnetic field (Fig. 8). The characteristic growth 
rates of this mode are nevertheless large compared to S2 dis- 
turbances. 
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Fig. 8 Same as in Fig. 7 but for S 1 modes. 



The models for the laminar solar tachocline are valid 
for microscopic diffusivities or sufficiently low turbulent 
diffusivities up to about 10 7 cm 2 /s (Kitchatinov & Riidiger 
2006). According to Eq. (8), this means growth rates up 
to a <~ 10~ 3 . The laminar tachocline models have toroidal 
field amplitudes of about 200 G or b ~ 0.002. The max- 
imum field in the models is attained at those depths where 
latitudinal differential rotation is reduced to about 10%. Fig- 
ures 7 and 8 show that the axisymmetric models for the lam- 
inar solar tachocline are stable to nonaxisymmetric distur- 
bances. 



3.2.3 Do stars other than the Sun possess tachoclines? 

Observations of Barnes et al. (2005) and theoretical models 
of Kitchatinov & Riidiger (1999) and Kiiker & Stix (2001) 
suggest that the absolute value AO of differential rotation 
is almost independent of the rotation rate for stars of given 
mass (but increases with mass). Figure 9 reformulates the 
stability map for the (most easily excited) Al modes in 
terms of the magnetic field strength and angular velocity 
of a star of solar mass assuming AO = aO = const and 
using Eqs. (2) and (3) to recover B$. 

The toroidal field amplitude in the laminar tachocline 
models is controlled by AO and fluid density only. There- 
fore, the limit B ~ 200 G of the solar tachocline model 
does also apply to the stars of same mass but different ro- 
tation rates. The star moves down on Fig. 9 along a vertical 
line of B ~ 200 G as the star spins down with age. The 
toroidal field can be only mildly unstable (<r < 10~ 3 ) for 
the laminar tachocline model to apply. All solar mass stars 
with rotation periods shorter than about two months should 
thus possess tachoclines after Fig. 9. 

Acknowledgements. This work was supported by the Deutsche 
Forschungsgemeinschaft and by the Russian Foundation for Ba- 
sic Research (project 05-02-16326). 



A' 




10 100 1000 10000 

B (Gauss) 



Fig. 9 Isolines of small growth rates of Al modes for various 
rotation rates and toroidal field amplitudes. The vertical dashed 
line shows the evolutionary track of a star according to the laminar 
tachocline models. Solar mass stars rotating not much slower than 
the Sun are predicted to possess tachoclines. 
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A Linear stability equations 

The equations used in our stability analysis represent the eigen- 
value problem resulting from the exponential time dependence 
exp(— iooi) of the disturbances. The dependencies on azimuth and 
radius are taken as Fourier modes exp(im4> + ikr). The equations 
are formulated in terms of scalar potentials for the vector magnetic 
and velocity fields (Kitchatinov & Riidiger 2007). 

The equation for the potential V of the poloidal flow reads 

a (cv) = -a 2 (is)-i^(iv) 



- 2^0, (CW) -2(1- n 2 ) 
+ 2 M fi A (CB) +2(l- M 2 ) 



2m 2 ^W 



d (fin) ow 

dfi dfi "" dfi 

d (h&a) dB 
dfi dfi 



9 (9fU - i- \ d(f£lA) 
+ 2m rl_2. B - mn A (CA) - 2m— ^ '-A 

dfi v ' dfi 

9 (f-ifi) , 2N dfi dV 



(Al) 



where A is the normalized radial wave length (5). The dimension- 
less eigenvalue u>, the angular velocity Q and the Alfven angular 
frequency CIa are all normalized to Qq, fi = cos 6, and 



The complete system of five equations also includes the equation 
for the potential W of toroidal flow 

w (CW) = -i ^ (CW) + mCl (CW) 
A 2 

d 2 

- mfl A (CB) - ((l - fi 2 ) fi) 

+ ™B^((i-S)n A ) + (iv)JL((i-» 2 )fi) 



(^)^((i-/* a )^) 



(A2) 



the equation for the potential B of toroidal magnetic disturbances 

cu (CB) = -i ^ (CB) + mi (ClB) - mC (Cl A W) 
A 2 

2dQ_ A _ d_ ( (l _ 2\2d&dA 
m dfi dfi I ^ ^ ' dfi dfi 



| m 2 ^A y — ( (i 2\2 d^A dV 
dft dft \^ 'dfi dfi 

the equation for the poloidal magnetic disturbances 

Co (CA) = -i (iA) + mQ (CA) - m£l A (CV) , 
A 2 

and the equation for the entropy perturbations 

ujS = -i %S + mClS + iV. 
A 2 

The equations include finite diffusion via the parameters 

vN 2 r)N 2 xN 2 



ngr 2 ' " 



ttgr 2 ' 



x ngr a ' 



(A3) 



(A4) 



(A5) 



(A6) 



where v, r\, and \ Me me viscosity, the magnetic diffusivity and 
the heat conductivity, respectively. 

The disturbances in physical units can be recovered with 



iCpiV 2 



P v = (n r 2 /k) V, T v = (Q r 2 ) W, a 



Pm = (v fiopQor 2 /k) A, T m = (y/fiopSlor 2 ) B, 



gkr 



(A7) 



and the vector magnetic and velocity fields can then be restored 
with 



r 2 r \ sm 



I i)T, u . d 2 P m \ 



+ 



drd6 ) 



r 



dT m 



1 d 2 P m \ 
sin 8 drd(f> J 



r 2 v r ysinO d(f> drd6 J 

+ e ± fdT 1 _ 1 d 2 P v \ 
r \ 80 sin 9 drd(f> J 

where e r , eg and are unit vectors in radial, meridional and 
azimuthal directions, respectively. 



(A8) 



T^-^d-p 



© 2007 WILEY- VCH Verlag GmbH & Co. KGaA, Weinhdm 



www.an-journal.org 



